LIQUIDATION IN LIMIT ORDER BOOKS WITH CONTROLLED INTENSITY 
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Abstract. We consider a framework for solving optimal liquidation problems in limit order books. 
In particular, order arrivals are modeled as a point process whose intensity depends on the liquida- 
tion price. We set up a stochastic control problem in which the goal is to maximize the expected 
revenue from liquidating the entire position held. We solve this optimal liquidation problem for 
power-law and exponential-decay order book models explicitly and discuss several extensions. We 
also consider the continuous selling (or fluid) limit when the trading units are ever smaller and the 
intensity is ever larger. This limit provides an analytical approximation to the value function and 
the optimal solution. Using techniques from viscosity solutions we show that the discrete state 
problem and its optimal solution converge to the corresponding quantities in the continuous selling 
limit uniformly on compacts. 

1. Introduction 

Liquidation of large securities positions has emerged as an important problem in financial math- 
ematics, linking together models of market microstructure and control theory. In this paper we 
consider an investor who liquidates a position through limit orders placed in a limit order book 
(LOB). The investor does so by choosing the price of the limit order; the higher the price of the 
limit order, the smaller the probability that it would be filled. The objective of the investor is to 
come up with an optimal limit order strategy that maximizes her expected revenue by date T. 

Our model for the above problem is based on a point-process view of limit order books which 
treats liquidation as a sequence of discrete events, i.e. order matches. More precisely, we assume 
that the investor effectively controls the frequency of her trades by choosing the spread s above 
the current bid price Pj. The trade intensity is controlled as A(s) and when a trade occurs, the 
investor generates a liquidation profit of s. Similar setups have been proposed in [5], [HI [15] and 
rely essentially on a queueing system representation of LOB's. 

A crucial modeling difference is whether execution takes place through market or limit orders. If 
investor trades via market orders, she necessarily encounters price impact through "eating away" 
a portion of the LOB. The precise price impact depends on the shape of the LOB, as well as its 
resilience. Conversely, there is no transactions or fill risk as market orders execute instantaneously. 
This point of view is taken in e.g., [231 [U [2]. On the other hand, if the investor trades through 
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limit orders, liquidation depends on being "lifted" by a sufficiently large market order, leading to 
substantial fill risk that again depends on the shape and depth of the LOB. The fill risk is related 
to the concept of virtual price impact [28] and is the focus of our model here. Related approaches 
to trading via limit orders can be found in [5l [T^ [20l I21j . Also, in our previous work [12] we 
considered the same LOB as here but with an uncontrolled trade intensity and temporary price 
impact from order size. 

Ideally, a fully specified model will reconcile the two approaches above, as well as consider the 
underlying risk preferences of the investor. This is especially important for dealing with simulta- 
neous trading on multiple exchanges, see the very recent preprints \23\ [22], as well as Section [5.51 
below. A full treatment of this problem will be the subject of a separate paper. 

Our starting point is a discrete-state problem for an investor holding n shares of an illiquid asset, 
n G N. Since practically speaking n is often large (on the order of hundreds of thousands), we also 
investigate the fluid limit of our setup. On a technical level, the fluid limit provides asymptotic 
results for the discrete-state problem (see Remark 13. 3p . which is the main focus of our paper. 

On a formal level, our control problem is equivalent to a controlled death process and is closely 
related to fluid approximations of some queueing problems. We refer to [8l|9l [lOl [181 [13 [26] and 
references therein for the most relevant strand of this rich literature. In contrast with the previous 
literature, which uses probabilistic arguments, we utilize viscosity techniques to show convergence 
(both of the value functions and the corresponding optimal controls) from the discrete- to the 
continuous-state problems. 

An investor holds n shares of an asset. Let {Pt)t>o be the bid price process for the underlying 
asset. Let r > be the risk-free rate. We assume that e~^*Pt is a martingale with respect to the 
optimization measure P on a filtered probability space {^,J-, (Gt))- This assumption is consistent 
with standard market microstructure models, see e.g. [2]. Let Af be the (controlled) intensity of 
order fill, and let > be the spread between the bid price and the limit order of the investor. 
Denote by Nt the ^-adapted counting process of order fills and Tk the corresponding arrival times. 



We assume that the investor has a deadline date T < +oo by which all trades must be completed. 
Remaining shares are liquidated at zero profit at T. 

To introduce the liquidation control, we assume that A^ = A(s(), so that the intensity of order 
fills is a function of the offered spread above the bid price. Moreover, we assume that the bid price 
P is unaffected by the limit orders created via (s^). Since e~'^^Pt is a martingale, the first term in 
(jl.ip is independent of Tj. Indeed, IE[X^"=i e'^'^^Pr^] = hPq and we may ignore P in the subsequent 
analysis. 



k 

Then Nt - £ A, ds is a martingale and expected revenue is 



(1.1) 
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(1.2) 



y(n,T):= sup E Y.""'"^^' 'nhn<T} 



(1.3) 



= sup E / e-^'^tclNt 



sup E 

ist)eST 



L 



■TAt{X) 



e-'^'stA{st)dt 



where 



t{X) := inf{t > Q : Xt = 0} 



is the time of liquidation. Here, Xt := Xq — Nt, with Xq = n, is a "death" (or inventory) process 
with intensity A{st). Note that T in (|1.2p represents time-to-maturity and St is the collection 
of J^-adapted controls, st > with J-'t := ct(A^s : s < t). The boundary conditions on V are 
V{n, 0) = Vn (terminal condition in time) and V{0, T) = V T (exhaustion). 

Remark 1.1. Our model is related to the limit order book setup of [5], which assumes that limit 
orders are "lifted" through sufficiently large market buy orders. Namely, a market buy order of 
size q, hits all limit sell orders that are within I{q) of the best bid. Assuming that buy market 
orders arrive in the form of a Poisson random measure on x with arrival intensity \dt 
and volume (mark) distribution f{dq), q > 0, a sell limit order at a given spread u is lifted with 
probability P(I(g) > u). By the thinning lemma on Poisson processes, such matching buy orders 
form a Poisson process with intensity A /^""i^^^ f{dq)- Empirical studies suggest a power-law depth 
function f{dq) oc q~^~°'dq [5] and therefore if A(0) < cxo, we can view our model within the [5] 
framework, A(s) oc [/"""^(s)]"", with the virtual price impact function |28j . 

As in [11], the above control problem can be transformed into a discrete-time Markov decision 
problem and the classical results from |13[ Ch. 8] can be used to prove a dynamic programming 
principle. Using the latter result one can show that the value function is a viscosity solution of 



with boundary conditions y(0, T) = F(n, 0) = and Vt denoting partial derivative with respect 
to time-to-expiration. Standard results also imply that an optimal control can be taken of Markov 
feedback type, = s{Xl,T — t). However, in most of the examples below we will obtain explicit 
solutions to this dynamic programming equation. Then a verification lemma can be used to justify 
that the solution of (|1.4p is indeed the value function. 

The optimization problem described in (jl.2p is simplified but highly tractable. In most of the 
examples below, we are able to obtain closed-form solutions which provide direct insight into the 
relationship between the LOB model and its depth function and the investor's liquidation strategy. 
In Section [2] we give an explicit solution for (jl.2p in the case of a power-law intensity control A(s). 
Section E] then studies convergence of the discrete problem (|1.4p to its continuous-state fluid limit. 
Our key Theorem 13.11 complemented by Proposition 13.21 and Corollarv 13.31 gives a full account of 
this convergence using techniques from viscosity solutions of nonlinear partial differential equations. 



(1.4) 



Vt + sup A(s) [V{n - 1, T) - V{n, T) + s] - rV{n, T) = 



s>0 
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In Section m we return to (jl.2p for the case where A(s) is of exponential shape; we are again able to 
provide several closed- form solutions. Finally, Section [5] considers several extensions and numerical 
illustrations of (jl.2p . including generic A(s), which shed additional light on the problem structure. 

2. Power-Law Limit Order Books 

In this section we assume that incoming buy orders have a power-law distribution for the spread, 
A(s) = for some q > 1. It can be observed from the computations below that if a < 1, then no 
optimal control exists. Similar assumption was made (and justified empirically) by [5] who write 
that in realistic markets a G [1.5,3]. 

Proposition 2.1. Assume that A(s) = As"" with boundary conditions V{0,T) = V{n,0) = for 
all n. Then the solution of (|1.4p and the optimal spread are respectively 

(2.1) y(n,r) = c„(l-e— ^)^/", .*(n,T)= f^y^^""'^(l-e— ^)^/", 

\arcnj 

with On satisfying the recursion 

(2.2) rcn = AaX{cn - c„_i)^"°, n > 1, Co = 0, 
where 

(2.3) A^ 



Remark 2.1. Note that V is "concave" in n in the sense that 

V{n+1,T) - V{n,T) < V{n,T) - V{n - l,r), 

i.e., its linear interpolation in n is concave in the usual sense. This follows immediately from (j2.2p 
since 

where the inequality follows from the fact that c„ (or V{n, T)) is increasing in n. The latter follows 
directly from (|1.2p . 

Also observe from (|2T]) and (12:2]) that 



s*(n,r) = ^{y{n,T) - V{n - 1,T)), 
a — 1 

which implies that n i— )■ s*{n,T) is a decreasing function, because n i— )• V{n,T) is "concave". One 
can also think of s* as the derivative of the linear interpolation of V in n. 

Proof. With the power law assumption ()1.4p reduces to 

(2.4) ~Vt + sup — (T/(n - 1, T) - V{n, T) + s) - rV = 0, 

s S° 
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and therefore the candidate optimal poUcy is 1 1— )• s*{Xt,T — t) in which s*(ri, T) = -^^^(y{n, T) — 
V{n — 1,T)). To begin solving this equation, we start with n = 1. Since y(0,T) = for all T, we 
obtain ioxV = V{l,T) 

-Vt + ^aAF^"° - = 0. 
This is a separable ordinary differential equation (ODE) which simplifies to 

T + C= — —dV. 



Using the boundary condition V{1, 0) = we integrate to obtain 

log{A^X - rV) = -raiT + C) ^ V{1,T) = !. ^{1 - e 



-raT^ 



Considering the equation for general n > 1 we therefore make the ansatz V{n^ T) = c„(l— e 

and plugging into (|2.4p the relation (|2.2p follows. □ 

Remark 2.2. If no discounting is present r = 0, one can verify that the solution of (|1.4|) is V{n, T) = 
dnT^^", where the sequence satisfies the recursion 

dn = X (^^) (dn - dn-l)'-'', do = 0. 



This result can also be obtained by taking the limit r — )• in (12. ip . ()2.2p . 

Fixing Xq = n, the inter-trade intervals cjj := Ti — rj_i, i < n have survival functions given by 

F{ai > t|r,_i) = exp ^-^ A(s*(n - i + 1, T - r,_i - s)) dsj . 

Noting that A{s*{n, T)) = ^^J^}^- for some constant C(n), it follows that A(s*(n, T)) dT = +oo 
for all n and e and therefore P(o"i < T — Tj_i) = 1 for all i < n. We conclude that even though there 
is no direct penalty if some orders remain at T, with probability one, the full inventory is liquidated 
by T, = P-a.s. In particular, the problem with a hard liquidation constraint V{x, 0) = 
— M1|^>0} for any liquidation penalty M > will have the same solution as in Proposition 12.11 

2.1. Infinite Horizon. As the execution horizon T grows, the investor faces a weaker liquidation 
constraint. Nevertheless, she still prefers to sell earlier than later due to the discount parameter 
r that incentivizes faster liquidation. For the limit T — )• oo we obtain an infinite-horizon model 
whereby strategies are time-homogenous. 

Taking T oo in we find that V{n) = Cn and s*{n) = A^/("-^)(arc„)^/(^-"). To under- 
stand how quickly execution takes place let us introduce expected time to liquidate S{n) which is 
defined to be S{n) := E[t(X*)|Xq = n], in which X* is the death process whose intensity at time 
t is A{s*{X^)) (representing optimally controlled inventory at t). When the inventory is X^ = n. 
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liquidation occurs at rate A(s*(n)), so that the interval until the next trade has an Exponential 
distribution with mean 1/A(s*(n)). It follows that 



(2.5) 



S{n) = Y^A{s*{j))-^ = A^/("-i) ^(Qrcj)-"/("-^). 



3. Continuous selling limit 

To better understand the results of Proposition 12.11 we consider a limiting continuous model. 
Let us denote the number of shares initially held by x. 

We first introduce a sequence of discrete control problems that converge to the continuous selling 
limit. For < A < 1, consider the problem where shares are sold at A increments and the intensity 
of order fills is A'^{s) := A(s)/A. We will denote by the "death" process with this intensity 
and decrements of size A. Then the resulting value function 



(3.1) V^{x,T) := sup E 



'x/A 

4 = 1 



sup E 

{st)€ST 



X G {0, A, 2A, • • • }, T G M+ would satisfy 

A 



(3.2) 



Vt + sup — :r{V^{x - A, T) - V^{x, T) + sA) - rV^ = 
s>o s"A 



in viscosity sense. 

Let us consider the first order partial differential equation (PDE) 



(3.3) 

which can be written as 



— Vx 

Vt + sup A rv = {), 



-VT + Ao,\vl " - rz; = 0, 



with boundary conditions f(x,0) = f(0,T) = 0. The solution of ()3.3p has the following determin- 
istic control representation 

^TAr(X(0).-) ^ 

{st)eST ■ 



(3.4) 



v{x,T) = sup / -^e~''^dt, 
{st)eST Jo St 



where dX- 



(0),x 



As J " dt, Xq^^'^ = X. Li fact, the solution of (j3.3p is explicitly given by 



(3.5) v[x,T) = (—\ ' x^"-^)/" (1 

\ra ) 

We denote the optimizer in (|3.3p by s('^)(x,T), which is explicitly given by 
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Remark 3.1. Plugging the optimizer back into the dynamics for X^^^'^ we obtain that 



dX-i ' = — Trdt^ 

* I — Q-ra(T-t) ' 



which can be explicitly solved as 



.(o),x I r ar 



(3.6) X, ' =xexp(^-y^ 

Let X £ M+ be fixed and let us consider all the collections {0, A, 2A, • • • } of grids that contain x 
as an element. In the next result, will show that as A — t- then V^{x) — t- v{x). In fact, the next 
result shows that this convergence is uniform on compacts. 

Theorem 3.1. As A ^ 0, — t- v uniformly on compact sets. 

Proof. Let us consider the regularized stochastic control problem 



(3.7) V^^''{x,T) := sup E 



X-/A 

E 

1=1 



e '^^'A-s,^l|,^<r} 



x G {0,A,2A, •••}, 



where := {s € St ■ st G [1/^, k]}, k > 1. Using a representation similar to the one in ()3.ip and 
using the lower bound on the controls s G S^, it can be seen that 

(3.8) V^'''(x,T) < -k'^-K 

r 

We will follow the arguments of [7j in the proof of their Theorem 2.1 (also see Theorem 4.1 
on page 334 of [19]) to show that V^'^ converges uniformly on compacts to the unique viscosity 
solution of 

k 

(3.9) -4+ sup A^^-TO*^ = 0, v''{x,0) = 0. 

s&[l/k,k] 

Let v'' and v'' be defined by: 

v''{x,T) := limsuplimsupsup|y^'*-'(y,5) : \x - y\ + \T - S\ < 6, y G {0, A, • • • }| , 

/(2;,T) := lim inf liminf inf fy^''=(y, 5) : \x - y\ + \T - S\ < 5, y G {0, A, • • • }| . 

By definition we have that vj' < < and that is lower semi-continuous, and is upper 
semi-continuous; see e.g. Proposition 5.2.1 of [6]. We will show that u*^ is a subsolution and that 
v'' is a supersolution of (j3.9p . It follows from Theorem 5.4.20 in [6] that a comparison result holds 
for this PDE (the compactness of the control space is required in order to apply this result). This 
comparison theorem would then imply that v'^ < . As a result, = = is the unique 
continuous viscosity solution of (j3.9p . This fact together with the way the functions and are 
defined also imply the local uniform convergence of V^'^ to . (For a similar argument see page 
35 of [17].) 

We now prove that is a viscosity subsolution of (j3.9p : the fact that is a viscosity su- 
persolution follows similarly. Let {xq^Tq) be a local maximum — (p for some test function 
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(j) G C^'^. Without loss of generality, we will assume that {xq,Tq) is a strict local maximum and 
that v^{xq,Tq) = <j){xQ^TQ), and 4) > 2^A;"~^ outside the ball B(xo,2o;r), where r > is chosen so 
that {xq,Tq) is the maximum of v'' — cf) on B{x(j,To;x). Thanks to the choice of the test function 
outside this ball, {xq,Tq) is in fact a global maximum of the function — cf) and it is attained on 
B{xq,Tq]x). (This is where the uniform boundedness assumption in (j3.8|) is used.) 

Let [x^ ,T^) G {0, A, • • • } X M_|. be a point at which V^'^ — 4> attains its (global) maximum. It 
follows from the definition of and the fact that {xq,Tq) is a strict global maximum of — (j) 
that there exists a sequence A„ — )• such that (x'^",T^") — )• {xq.Tq), V^"'^ — (p attains its global 
maximum at that point and y^"''=(x^" , T^") — )• v^{xq,Tq). From the global maximality 

y^"'''(x,r) - yA„,fc(^A„^^A„^ < (t){x,T) - (l){x^",T^"). 

Moreover, it can be argued as in [TT] using the discrete dynamic programming principle (see |13j ) 
that V'^"''' satisfies 

-V^"''' + sup (v^-'''{x^- - An,T^-) - V^"^^{x^-,T^-) + sA„) - rV^"^'' = 

se[i/fc,fc] s^An V / 

in the viscosity sense. Then 

-4)T+ sup ((/)(x^" -A„,r^")-(/.(2;^",r^") + sA„) -r(/> + r((/>-y^"''=) >0. 

Taking the limit as A„ — )• we obtain from this equation that 

-(j)T{xQ,TQ) + sup ^{s-(l)x{xQ,TQ))-r(t){xQ,TQ)>^, 

sg[l/fc,fc] S 

which proves the subsolution property of . Here, we exchange the limit in A„ and the supremum 
with respect to s using Proposition 7.32 in [13] which we can apply thanks to the compactness of 
the control space. 

It follows again from Theorem 5.4.20 in [6] that the unique solution of ()3.9p is given by 

(3.10) v''{x,T)= sup / 7- — - — - — -dt, 

^ ' ^ ' is,)eSrJo ((st V 1/fc) A fe)"-i 

where dX^^^'^ = —\/{{st V 1/k) A k)°'dt, X^^^'^ = x. We will show that converges pointwise to 
v: 



hm V (x, T) = sup sup / —, — — — — ^ dt 

k is,)isrJo ((st V l/fc) A /c)"-l 



TAr(X('=)'^) Ae"^'* rThr{X) ^ 



-rt 



= sup sup/ w 1 // \ A u\a~i - / -^dt = v{x,T), 

{st)eST k Jo ((st V 1/A:) A /c)" (s^e^T ^ 4 

where the inequality follows from the lower semi-continuity of the map X i— ?• t{X); see Lemma 5 in 
|18j . On the other hand, since u is a supersolution of ()3.9p . the comparison result Theorem 5.4.20 
in [6] implies that <v for each k, and as a result 

lim v^{x,T) < v{x,T). 

k-^oo 
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Combining the last two inequalities, we obtain the pointwise convergence of v'' to v. Pointwise 
convergence, on the other hand, implies uniform convergence on compacts due to Dini's theorem, 
since we already know that f is a continuous function of its arguments, and that v'' is an increasing 
sequence of functions. The latter fact follows from the fact that is a supersolution of the PDE 
v'' satisfies. □ 

Remark 3.2. Results somewhat similar to Theorem 13.11 appeared in [Hj \9\ IIOI 118^ 1251 [26] which are 
on the optimal control of queueing networks. (Among these papers only [18] considered optimal 
time-to-empty queueing control problems.) To prove Theorem 13.11 we used a completely different 
approach than the above literature, which had relied on probabilistic arguments. Our approach 
relies in contrast on the analytical approximation ideas of We see the prelimit control problem 
as the discretization (only in the space variable but not in the time variable) of the "fluid limit" 
first order non-linear PDE (j3.3p and rely on convergence of the approximation schemes to the 
viscosity solutions of such non-linear PDEs. This approach could be fruitful in general in proving 
"fluid limit" results associated to controlled queueing networks. 

The following is a strengthening of Theorem 13.11 which is an interesting result in its own right. 

Proposition 3.2. For any sequence (A^) with = 62~^ , we have V^'^ '\ v as — )■ oo. 

Proof. We show that for any A > 0, V'^'^ < V^. Due to the factoring of T and x in Proposition 
12.11 it suffices to establish this result on the infinite horizon where strategies are constant between 
trading times. 

Fix e > and let s^'^ be an e-optimal strategy for V"^^. This policy is defined over x G 
{0, 2A, 4A, • • • }. We will recursively construct a policy over the domain x G {0, A, 2A, • • • } 
that outperforms s^^. The dynamic programming principle implies that 



V^^{2nA) < E e-^^i [2s^^(2nA)A + y^^((2n - 2)A 



rn [o ^2 A (-9^ A \ A , T a2A / 



+ e 



Similarly, given the liquidation strategy s'^ and corresponding trading times fj, the resulting ex- 
pected profits denoted as V^{x), x E {0, A, 2A, • • • } satisfy for y = 2nA, 



V^{y) = E [e~^^is^(y)A + e-^^^ j^A^^ _ ^ yA^y _ 2^)| 



2A(s^(y))+r\ 2A(s^(y - A)) + r 

Given S2n = s^^(2nA ) we prove below that there exists n G M+ , such that 

A(s2„) , , , 2A(n) r , 2A(n) , , 
3.11 w (2s2nA + V)< , , , <^ ^A + ' iuA + V) 

^ ' A(s2n) + r^ ^-2A(n) + r\ 2A(n) + r^ ' 

for any F > 0. This would establish y^(2nA) > y'^(2nA) > V'^^{2nA)-£ by induction on n after 
setting s'^(2nA) = s'^((2n-l)A) = u. Since e is arbitrary, the statement of the proposition would 
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then follow. Note that in the above construction, the A-investor trading in smaller increments and 
twice as much, uses the same spread u to trade when her inventory is 2nA or (2n — 1)A. 



Let := j^^^^'-^^^ and define u implicitly through 2A{u)+r '■= ^ < ^- Solving for S2n and u in 
terms of z and using A(s) = As^" we obtain 



y rz^ J \ rz 

Observe that by construction 

A(s2n) 4A(n)2 



A(s2„) + r (2A(«) + r)2' 
so that the Laplace transform at r of the duration to execute two trades by the A-investor is equal 
to the Laplace transform at r of the duration to execute one trade by the 2A-investor. Using this 
fact, (|3.1ip is equivalent to 

2A(s2n) ^ 2A(n) / 2A(n) 

-S2n < ^ w N i + 



A(s2n) + r ~ 2A(m) + r { 2A{u) + r 

2z'^ S2n < z{l + z)u. 

Since a > 1 and the terms on both sides of the above inequality are positive, we may raise both 
sides to the a-power and plug-in the expressions for S2n and u to find 

rz^ rz 
= 2\r-\l - z){l + z)z"-H(2z)°-^ - (1 + z)"-i] < 0, 

where the last inequality follows since z < 1 and a > 1. This shows that (j3.1ip holds and concludes 
the proof of the proposition. □ 

Next, we show that the strategies also converge thanks to the concavity of all the functions 
involved. 

Corollary 3.3. Let us denote by s^^'^ the pointwise optimizer in ()3.2p . Then we have that 
s^^\x,T) —7- s^^\x,T). (Here, x is fixed and we take the limit over the grids that pass through x.) 

Proof. On the one hand, as in Remark 12.11 we can think of s^^^ as something proportional to the 
left derivative of the linear interpolation of , denoted by V^, which is increasing and concave. 
On the other hand, s^^'^ is proportional to the derivative of the concave difFerentiable function v. 
By Theorem 24.5 on page 233 of [27] it however follows that for any a; > 0, 

\D-V^{x,T)-v^{x,T)\ <e 

for small enough A, where denotes the left derivative operator with respect to x. □ 

Remark 3.3. Theorem 13. II tells us about the asymptotics of c„ in (|2.2p : 

A \ 

— 1 as n ^ oo. 

ra 
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Corollary 13.31 can be used to find out the marginal price asymptotics in Proposition 12.11 

(3.12) s*(n,T) ~ ( — I —rr- as n — )■ oo. 

^ ' ^ ^ \ar J nV" 

Clearly, the spread will go to zero as n — t- oo. But here we are able to obtain the rate of convergence 
to zero as a function of the remaining inventory. 

For time to execution on infinite horizon we have for ri = inf{t : Xt < xi} and S{xi,X2) '■= 
E[ti|Xo = 2:2] that 

(3.13) S{xi,X2)= , . du, 

since intuitively when inventory is of size u, the expected time to liquidate an infinitesimal quantity 
du is inversely proportional to the current trading rate A(s('^) (^))- Plugging in s^^^ (n) = {-^^)^^"' we 
obtain A{s^^\u)) = aru or S'(xi,X2) = ^log(^) which shows that orders are filled in logarithmic 
time (as xi ^ the remainder is executed arbitrarily slow). 

4. Exponential-Decay Order Books 

The power-law order book implies that trades can be made arbitrarily quickly as the spread goes 
to zero: lims_^o A{s) = +00. Also, it gives a relatively good chance of executing trades deep in the 
book, i.e. when s is large. For less liquid markets, both of these features might not be realistic. 
Accordingly, we consider an exponential-decay LOB, with 

(4.1) A(s) = \e-^', K > 0, 

where k controls the exponential depth of the book and A = A(0) is the order intensity at the bid 
price. The optimization problem for the spread is now of the form 

sup Ae"^" {V{n - 1) - V{n) + s) , 

s>0 

which leads to the candidate optimizer s*(n) = + (^(^) ~ V{n— 1)). We observe that s* is 
bounded away from zero so no trades are ever placed close to the bid. 

4.1. Finite Horizon. With a finite horizon and no discounting we obtain the following closed- form 
solutions to the execution problem. 

Proposition 4.1. Consider again V^{x,T) defined in (j3.ip with boundary condition V^{x,0) = 
V^{0,T) = 0, r = and A(s) given in ()4.1|] . Then for x = nA, 

(4.2) K^(,r,.^,o.(tK0)' 

and 

.*(nA,T) = l 1 + log 1+ 



^j=0 (Ae)J j! ^ 
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As A \ 0, V^{nA,T) — )■ v{x,T) uniformly on compacts, where v{x,T) solves the nonlinear first 
order PDE 

A 



(4.3) VT{x,T) = -e 



-1-KV^(X,T) 



with boundary conditions v{0,T) = v{x,0) = 0. The solution to this PDE satisfies 

(4.4) - log (-t] < v{x, T) < —T. 

K \X J KC 

Remark 4.1. In the above notation x is the number of shares, which is fixed across the problems. 
When we are taking the continuous liquidation limit, we let A J, 0, the size of trading units, while 
taking the number of units n as nA = x, for x constant. 

Proof. V^{nA,T) satisfies the HJB equation 

-drV^iriA, T) + sup XA-^e-'^'iV^iin - 1)A, T) - V^{n, T) + sA) = 0, 

s>0 

which can be written as 

(4.5) dTV^{nA,T) = ^ exp (-1 + A~i/t(y^((n - 1)A, T) - y^(n, T))) . 
Letting B := integrating and using V^{-, 0) = we find for n = 1 that 

(4.6) V^{A,T) = —log{l + BkA-^T). 

K 

Iterating over n the separable ODE for V'^{n,-) in (|4.5p we obtain (|4.2|) . The expression for the 
optimal spread follows from s*{nA,T) = ^ + ^ ((n-i)A) ^ 

The proof that as A J, 0, V'^{nA,T) — ?• v{x,T) uniformly on compacts can be proven as in 
Theorem 13.11 

Observe that both bounds in (14. 4p satisfy (j4.3p . To prove the lower bound, let us introduce the 
following function: 

Clearly, < V^. From Stirling's formula we know that 

where we use ~ to indicate that the ratio of the left to the right-hand-side converges to 1 as 
n —7- oo. As a result, y^(a;,T) ~ x / K\og{\T / x) , recalling that x = nA. Now the lower bound 
in (j4.4p follows since v{x,T) > limA4,o^'^- We could have provided an alternative proof using 
a comparison theorem for the first-order non-linear PDE (j4.3p since in fact x/Klog{XT/x) also 
satisfies this PDE with a smaller boundary value at T = 0. We preferred to be more constructive 
in our proof. 

The fact that XT/(Ke) is an upper bound on v follows directly from the observation that vt < 
X/{ne) (recalling that v is increasing in x) and that v{x,0) = 0. 

□ 
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Remark 4.2. Since the trading rate is bounded A(s*) < Ae^"*^, for x > Xe^^T the fuh inventory 
cannot be hquidated by horizon T. Therefore, in the region V := {{x,T) : x > Ae~^T}, v is 
independent of x and the upper bound is tight: v{x,T) = -^T on P. 

Remark 4.3. Here, we will determine the shape of t — t- xj:^^'^ in Remark 13.11 for the exponential 
order books. First, 

(4.7) dxf ^''^ = -A(s(xf T - t))dt, 

where A is given by dH]). On the other hand, s(xf T - t) = i + 't;^(xf T - t). Using this 
relationship, along with (|4.3p . which implies that 1 1— )• s{X^^^'^ ,T — t) is a constant function (let us 
denote that value by s*), it follows from ()4.7p that 



0, tG[0,T], 



i.e., t ^ Xf^'"", t £ [0,r], is a strictly decreasing linear function. In fact, one can compute s* by 
maximizing the value function (j3.4p (after replacing power rate with exponential) over constant 
spreads (since the optimal spread is known to be a constant). This yields that s* = log (^) if 
XT/x > e. Otherwise s* = The expression for the optimal spread and (|4.7p in turn imply that 
if X < Ae-^T, = x(l - ^) (which is at T) and if x > Xe~^T we have that ^'"^ = x- Xt/e 

(which remains strictly positive at T). 

4.2. Infinite Horizon. We also have closed-form expressions for the infinite horizon case. 
Proposition 4.2. For exponential- decay LOB with T = +oo and discounting rate r > we have 

(4.8) V^{x) = (^Ar-^A-^exp(/^ ^^^^~^^ - l)^ , x £ {0,A,---}, V^{0) = 0, 

where W is the Lambert- W function (or the double-log function defined as z = W(y) for ze^ = y)- 
As A ^ 0, V^{x) — )• v{x) uniformly on compacts where 

(4.9) ^^^enrv{x)^^_e_rx^ 

and li{y) := j^dt is the logarithmic integral function. 
Proof. The HJB equation for y^(a;) is 

-rV^{x) + sup AA-^e"'^" {V^{x - A) - V^ix) + sA) = 0. 

Using the optimizer s*^'^^ = « + ^ (^~^) -^^g reduce to 

, _ AA f ^ , _y^(a;- A) -"l^^(a;] 



V (x) = exp [ —1 + K 



Kr \ A 
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which has closed-form solution given by (j4.8p . Arguments similar to Theorem 13.11 imply that 
— ?• V uniformly on compacts and the continuous inventory limit satisfies 

(4.10) -rv{x) + sup Xe^'^'is - v'{x)) = 0. 

s>0 

Solving for v' we obtain 

\ If-, 1 / nrvix) \ '] 
„'(x) = --|l + log(^j|. 

The last nonlinear first-order ordinary differential equation has closed- form solution given in ()4.9p . 
Asymptotically Imix^^ v{x) = and the optimal spread is 

□ 

We note that since v is increasing in x, x i— t- s^^\x) in ()4.1ip is decreasing and so v is concave. 
As before, \\.m.x^Q s^^\x) = +oo so the control space remains unbounded, however the pay-off rate 
s^'^^A(s^*^^) is bounded. Moreover, a direct check verifies that and v are inversely proportional 
to the exponential depth parameter k, i.e. doubling k (making the order book more shallow) halves 
and u, and correspondingly halves the optimal spreads s^^'^ and s^'^^ 

Figure [1] graphically illustrates the difference between exponential-decay and power-law LOB's. 
As observed, for an exponential LOB, s* is bounded away from zero, while \\mx^ooS*{x) = in 
power LOB's. Moreover, while YiTa.x\.o s^^\x) = +oo in any LOB, the rate is much slower in an 
exponential LOB (due to thinner tail for large spreads) compared to power-LOB. 



5. Discussion and Further Extensions 

5.1. Numerical Example: Convergence to the Fluid Limit. To illustrate the convergence 
to the fluid limit consider the problem of selling up to rr = 5 blocks of shares in a power-law 
LOB with A(s) = s"^. We suppose that a block corresponds to 100 shares and that the minimal 
trading unit is either 5 or 1 shares, i.e. A = 0.05 and A = 0.01 respectively. For a fixed A, we can 
easily compute V^{x = nA) or v{x) using the results in Section[2l Figure [2] illustrates the percent 
difference between and the fluid limit v. As shown in Proposition 13.21 is decreasing in A 
and limAj^o = We observe that the convergence is quite rapid in x. 

The right panel of Figure [2] shows that the controls themselves are also very close. We observe 
that s^^'^ \ s^^\ Note that here we are essentially comparing s^*^^ with its right-sided Riemann- 
sum approximation, since s^^\x) corresponds to the spread charged for all shares in [x,x — A) 
while s^^\x) corresponds to the marginal spread at x. Accordingly, better approximations, such 
as s^^\x) := j^-A •s^'^H^) du, would make the discrete and fluid controls even closer. 
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Figure 1. Optimal controls for power-law and exponential-decay order books. We 
take r = 0.1 and depth functions A(s) G {s~'^,s~'^,e^~^}, which have been normal- 
ized such that A(l) = 1 in all three cases. The plot shows the resulting fluid limit 
spreads s^^\x). 





v(x)/V'^-'''{x) 

V0.0^j.)/V0.05 ^^.J 


\ ^ 






y 




1 

1 




_ 1 
I 





0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 ' 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 



Figure 2. Convergence to the fluid limit. Left panel: the ratio between discrete 
and continuous V^{x)/v{x) for A = 0.05 and A = 0.01. Additionally, we plot 
yA^yA deflned in ()5.ip . Right panel: ratio of the fluid limit optimal control 
to the discrete s^^\x) for A G {0.01,0.05}. 



Given the simple expression for the fluid limit control s^^\x), a useful approximation is to use 
a discretized version of s^^^ as an approximately optimal control for V^. Let 



(5.1) 



V^{x) := E 



'x/A 

^e-'^"»sW(x-iA) 

i=l 



X G {0,A,2A,...}, 
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represent the expected gains from a discrete strategy which uses a spread of s^^\{n — i)A) for the 
i-th trade of size A. In the left panel of Figure [2] we see that this approximation is excellent for 
even for moderate values of x (less than 1% difference for A = 0.01 and x > 1). 

5.2. Execution Curves. A popular way of describing a trade execution algorithm is through 
the execution curve t i— )■ Xi/Xq, see e.g. [3l HI [20]. In our model with execution risk, Xt is a 
random variable, and we will therefore consider the natural analogue of average execution curve 
E{x,t) := E[Xj*'^], where X*'^ is the remaining inventory at t < T starting with initial condition 
Xq^ = x . The baseline case where E{x,t) = x(l — t/T) is linear, corresponds to "linear price 
impact" or zero-risk-aversion in [3] and implies that the average trading rate is constant. 
For notational convenience we temporarily fix A = 1. Recall that 

dX;'"" = -A{s*iX;''',T - t)) dt + dMt 

where (Mt) is a martingale (the compensated order departure process), which implies by an appli- 
cation of Ito's formula that {E{x,t))'^Q satisfies the system of inhomogenous linear ODEs 

^^^^ = A(s*(x,T - t)){E{x - - E{x,t)), E{x,0) = x, 

with E{0, t) = 0. Any finite collection of these ODEs can be solved analytically using integrating 
factors, or numerically with any standard solver. 




D 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0,1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Time ( Time ( 



Figure 3. Execution curves for different power-law books. We take Xq = 6, r = 
0.1, T = 1, and depth functions A(s) G {s"^,s~^}. Left panel: average inventory 
E{x, t) as a function of time t. Right panel: average trading rate as a function of 
time. 

Thus, whenever we have an explicit formula for the optimal spread s*{n,T), E(x,t) is also 
available analytically (or more practically as a solution of an ODE). For the case of power-law 
order books, we obtain from Proposition 12. II that 

A~^(Qrcfc)°/("-^) 



A{s*{k,T-t)) 



1 — exp(— ar(r — t)) 
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We observe that depending on the parameter values (in particular initial inventory k vis-a-vis order 
shape a), A(s*(A;,T)) could be smaller or bigger than k/T, i.e. the ordering between the initial 
trading rate and constant trading is ambiguous. At the other end, as i — >■ T, the spread s*{k,t) 
goes to zero and consequently limt^T dE{x,t)/dt = — oo. 

Figure [3] shows that as a increases, the shape of i i— )• E{x,t) changes substantially. In particular 
for large a (corresponding to "thinner" power laws), the execution curve has an S-shape, with 
trading rate high in the beginning and end of the time interval. On the other hand, for a small, 
the execution curve lies entirely above the baseline, i.e. the limit order trader consistently executes 
slower. This occurs due to the two competing effects of trying to extract profit (which slows 
execution) and the time decay, i.e. the need to make the deadline which speeds up trading. We 
observe that when the limit order book is deep (small a), the profit effect dominates; this is a new 
phenomenon compared to most existing models, such as |20j . 

Finally, as a comparison. Figure [3] also shows the deterministic case t i— )■ X^^^'^, see (j3.6p . In 
the latter case, 1 1— ?• xj.^^''^ is strictly convex, which contrasts strongly with the pre-limit situation. 
As A —7- 0, execution risk vanishes and the discounting effect takes over, making the investor sell 
more in the beginning. Indeed, in the fluid limit the investor can smoothly drive X^^^'^ to zero, 
while for A > 0, F{X^''-^ > 0) > for any t < T, but X^'"" = since dE{x,t)/dt\t=T- = -oo. 

In the exponential order book case with no discounting, we observe that the trading rate 
A(s*(n, T — t)) is independent of the depth parameter k. Moreover, numerical experiments suggest 
that E[x, t) is (slightly) convex in t, i.e. the trading rate is monotonically decreasing in time. This 
agrees with the classical results of [3] . Finally, we recall that Remark 14.31 shows that in the fluid 
limit, the execution rate is constant over time, and t i— )■ X^^^'^ is linear. This occurs because 
{X^^'^ ^t) is the characteristic curve of the PDE given by (j4.3p . This phenomenon resembles the 
constant trading rate in Alfonsi et al. [1] who studied (continuous) trading through market orders 
only, with the LOB depth function driving the price impact mechanism. Again, we find a sharp 
dichotomy between the deterministic limit where Xj^, = for x < Ae~ T and the stochastic 
version where E(x,T) > strictly for all x > 0. 

5.3. General Order Book Depth Functions. Our basic setting can be readily extended to 
allow for more sophisticated or complex models. Below we review several such extensions; for ease 
of presentation we treat them in the stationary infinite-horizon setting. 

Let us revisit the optimal execution problem for a generic order book depth function A(s). In 
general, there are no closed-form expressions for V{n) and the continuous fluid limit v{x) becomes 
a useful analytic tool to understand the solution structure. In that regard, both Theorem 13.11 and 
Corollary 13.31 continue to hold under some reasonable assumptions on the intensity function A. 

Theorem 5.1. Consider the optimal liquidation problem on infinite horizon with a general intensity 
depth function A. Then the statements of Theorem \3.1\ hold, and if we further assume that the 
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function x i— )■ A{x) is decreasing and that 

, , A(a;)A"(x) 
(5.2) (V(x))V 

then both and v are concave, the corresponding controls s*-'^^ and s^*^^ are decreasing and the 
conclusion of Corollary 13.31 still holds. 



Remark 5.1. For power-law LOB's, condition (j5.2p holds precisely when a > 1, while for exponen- 
tial LOB's it always holds. Both of these order books have decreasing intensity functions. 

Proof. The proof of Theorem 13.11 can be done without much change since we did not make use 
of any special properties of A(s) there. We will prove the stated concavity and monotonicity 
properties from which the statement of Corollarv 13.31 follows immediately as before. 

By time-stationarity between trading dates the controls are constant and the dynamic program- 
ming principle until the first jump time for V^{x) gives 

V^{x) = sup r Aifl e-(^(^)^"'+^)*(sA + V^{x - A)) dt = sup ^ , (sA + V^{x - A)) . 

s>o Jo A ^>o A(s) + r A ' 

Differentiating the right-hand-side with respect to s, the first order condition for s* = .s^^\x) is 
rA'(s*)(s*A + V^{x - A)) + A(s*)(A(s*) + rA) = 

^ ry^(x - A) = -rs*A - |;^(A(s*) + rA) := F{s*). 

is non-decreasing; therefore, if the derivative of F is negative, then s^^\x) decreases in x. 
Explicitly, 



(5.3) F'(s*) = -(rA + A(s*)) 



AA 



(A 



k{s*)k"{s*^ 
(A'(s*))2 



< ^ 2 > 



Thus, ()5.2p is sufficient for x i— )• s*{x) to be decreasing. Under this assumption and the assumption 

h(s*{x))+r 



that A is decreasing we would have that x i— )■ A(^*fa)y+rA increasing, and as a result 



V\x) - V\x - A) = + - A)) 



~ {s*{x - A)A + V^{x - 2A)) 



K{s*{x - A)) + rA 

< . {s*{x - A)A + V^{x - A) - s*(x - A)A - V^{x - 2A)} 

< V^{x- A)-V'^ix-2A), 
so that is concave. 

Concavity of on the other hand implies the concavity of v. This is thanks to the first assertion 
of the theorem from which we know that converges to v uniformly on compacts. Next we will 
show that s is a decreasing function. The value function v satisfies the first order PDE: 

supA{s){s-v')-rv = 0, ■u(O) = 0. 

s>0 
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Optimizing over s yields 

Hence v is concave if and only if the right-hand-side above is a decreasing function of x. However, 
it follows from (15.21) that the function 



My) 

is increasing. As a result the concavity of f , which we have already shown, is equivalent to 
X I—)- s^'^\x) decreasing. □ 

5.4. Regime Switching Market Liquidity. Empirical evidence suggests that market liquidity 
is not constant, |14] . As a first step towards capturing more complex liquidity behavior, we consider 
a simple regime-switching model for market activity level in the power-law LOB's. More precisely, 
suppose that arrival rates are modulated by a two-state Markov chain M with states {0, 1}, 
representing an active market and 1 representing a slow market. We will denote the transition rate 
from to 1 by and the transition rate from 1 to by ^i. Under regime 0, the arrival rates of the 
orders are Ao(s) = Aq/s", and under regime 1, the arrival rates of the orders are Ai(s) = Ai/s". 
We will take Aq > Ai. We assume that M is observed and known by market participants. 

Denote by U{n) (respectively W{n)) the infinite-horizon value function for an inventory of n 
shares under the active (resp. slow) market regime. The value functions satisfy the following system 
of equations: 

Ac,\o[U{n) - U{n - 1)]^-" - rU{n) + 0o[VF(n) - U{n)] = 0, 
AaXi[W{n) - W{7i - 1)]^-" - rW{n) + 0i[C/(n) - Win)] = 0, 

with terminal condition U (0) = W^(0) = 0. The continuous selling approximation of these functions, 
which we denote by u and w respectively satisfy the following system of ordinary differential 
equations: 

AaXoulr"' - (r + eo)u{x) + Oowix) = 0, 
AaXiwl'"" - (r + 9i)wix) + 6iu{x) = 0, 



(5.4) 



with u(0) = u;(0) = 0. 

Proposition 5.2. The solutions to (j5.4p are u{x) = c^x'^ and w{x) = c*x^ where p = (a — l)/a 
and 

(5.5) (^Y\cl<4<(^^'''^ 

\ra J \ra 

Observe that the bounds in ()5.5p correspond to the single-regime solutions given in (|3.5p with 
T = +00. 
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Proof. We begin with an ansatz of u{x) = cqx^ and w{x) = cix^ with p given in the statement of 
the proposition. Comparing with ()5.4p . the coefficients cq and ci need to satisfy 

A«AoP^-"4-° - + ^o)co + Ooci = 0, 

A„Aipi-"c}-° - (r + ^i)ci + 9iCo = 0. 

Re-writing as 

r + 6i Ai 

(56) = — ■ 

it easily follows that this system of equations has a unique solution (cg, c^). Indeed, cq as a function 
of ci is strictly increasing and goes from — cxo at ci = to oo and ci = oo. Similarly, ci as a function 
of Co is also strictly increasing from — oo at cq = to oo at cq = oo. It directly follows from these 
facts that these two curves intersect and do so at only one point. Moreover, the identity function 
Co = ci, intersects the first function (j5.6p first, and the second function in the same equation last. 
This proves the ordering in (|5.5p and concludes the proof. □ 



Proposition 15.21 shows that the asking spread will always be higher under the active market 
regime when the order book is deeper. 

Figure H] illustrates the impact of multiple liquidity regimes. We take Aq = 1.5, Ai = 0.5 so 
that trade intensity is tripled in the active regime. We plot c* as a function of ^o = = ^ for 
r = 0.1 and a = 2. As 6* 0, we have c* VoJX^, while as oo, c* the 
fast-switching limit. 




Figure 4. Regime switching model. We take oq = 2, Aq = 1.5, Ai = 0.5 and 
r = 0.1. The regime-switching rates are equal ^o = = ^- 
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5.5. Two-Exchange Multi-scale Model. Another possibiUty is to consider an investor trading 
on multiple venues. For example, suppose the investor can liquidate her holdings through two 
different exchanges, with each exchange possessing its own LOB. To distinguish the two exchanges, 
we suppose that on exchange C(ontinuous) the orders are infinitesimally small, but on exchange 
L(arge) they are of large but finite size relative to the total order size. More precisely, we assume 
that in the continuous limit, the exchange C orders are infinitesimal, but exchange L orders are of 
size 6. If the remaining inventory is less than x we assume that the next trade on exchange L will 
liquidate the entire x. In other words, actual trades on exchange L are of size mm{6,x). 

To keep the model tractable, we assume that each exchange has power-law depth with identical 
depth parameter a > 1. The resulting time-stationary value function v{x) solves 

(5.7) sup Ao h sup Ai rv = 0, v[0) = 0. 

so>0 Sq si>0 

Plugging in the first-order optimizers leads to 

(5.8) AaXov'{x)^-'' + AaXi{xA6)''{v{x)-v{{x-6)+))^-''-rv = 0, v{0) = 0. 
Lemma 5.3. There is a unique solution to (j5.8p . 



Proof. Equation (j5.8p is a first order nonlinear delay ODE and can be solved by successive patching. 
Namely, first solve the ODE 

(5.9) ^aAot;(o)(a;)^"" + AaAix"(i;(o)(x))i-" - rv^o) = 0, 
with W(o)(0) = on [0, 6]. We then solve 

(5.10) A«Ao^;[i)(2;)i-" + A«Ai5°(^;(i)(x) - v^o){x - -J))'"" - rv^i) = 

on [6, 26] with initial condition V(i)((^) = ^'(o)('^)- In (|5.10p we treat z;(o) as a source term, observing 
that f (0) {x — 6) with x G [5, 26] has already been computed before. Proceeding in this fashion, 
we finally set v{x) = t'(„)(x) for x E [n6, (n + 1)6] to recover the global solution. On each of the 
intervals [n5, {n + l)6] the corresponding ODE has a locally Lipschitz driver so classical results give 
existence/uniqueness of solution □ 

Remark 5.2. Numerical computation of the solution of (j5.9p should be handled with care since 
v|q^(0) = oo. We get around this singularity using the following observation: For x small enough, 
the benefit of large orders is negligible since the probability of getting a large order is very small. 
Therefore, close to zero, v{x) ~ ^0(3;) = from (|3.5p . 

We also remark that the solution to (j5.9p is in general no longer concave, with concavity likely 
to fail around the knots 6,26, .. ., where the derivative v' does not exist. 



Typically, trading intensity on the small-order exchange is several magnitudes larger than via 
the big trades (done through e.g. a proprietary dark pool, see e.g., |22j where a single large dark 
pool trade liquidates the entire position), so Aq ^ Ai. Fixing the time-scale as Aq = 1, we are 
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therefore led to consider an asymptotic expansion in small Ai. Formally, let Ai = Xe for e small 
and consider a power series expansion in e, 

v{x) = vo{x) + evi{x) + e'^V2{x) + . . . , 

Plugging into (j5.8p and matching powers of e we find that vq{x) solves the 1-exchange problem of 
(|33|) . so that ^0(2;) = Next, 



Aaa{l — a)rxv'i{x) + AaX{5 A xY 



^] " (x^ - (x - 6)7^y-'' - rv,{x) = 0. 



This is a first-order linear ODE with non-constant coefficients and therefore vi{x) can be expressed 
in closed-form using integrating factors as 

„2-a-l 

(5.11) vi{x) 



-Be 



C2y 



y 



X < 6, 
X > 6, 



with C 



Xar f Xo_\(l— a) / a _ 
(1— a) V or / 



Bc,-2+a 



-T, C2 = AaXS^'Ba 



{l-a)/a 



and Ba 



-. The latter 



integral only involves powers of y and can be easily computed numerically. Similarly, the equations 
for higher-order terms are again first-order linear ODEs and so V2, etc., can be written iteratively 
in closed-form. 
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